#=
相互作用
=#

using SpecialFunctions


"""
库仑势能
"""
function bare_Coulomb(r1, r2)
    dis = @. (r1 - r2)^2
    return 1 / sqrt(sum(dis))
end


"""
Kelbg势能
"""
function _Kelbg_potential(rij::Vector{T}, λab::Float64) where {T}
    absrij = sqrt(sum(rij.^2))
    rrel = absrij/λab
    c2 = -exp(-rrel^2)
    c3 = sqrt(π)*rrel*(1-erf(rrel))
    if isapprox(0.0, (1 + c2), atol=eps())
        return sqrt(π)/λab
    else
        return (1 + c2 + c3) / absrij
    end
end


"""
电子之间相互作用
"""
function Kelbg_repulse(r1::Vector{T}, r2::Vector{T}; λ11=_WLElec) where {T}
    return _Kelbg_potential(r1-r2, λ11)
    #dis = @. (r1 - r2)^2
    #rabs = sqrt(sum(dis))
    #rrel = rabs/λ11
    #c2 = -exp(-rrel^2)
    #c3 = sqrt(π)*rrel*(1-erf(rrel))
    #if isapprox(0.0, (1 + c2), atol=eps())
    #    return sqrt(π)/λ11
    #else
    #    return (1 + c2 + c3) / rabs
    #end
end


"""
所有粒子之间的相互作用
"""
function Kelbg_repulse_pairs(rpos::Matrix{T}; λ11=_WLElec) where {T}
    ϕ = 0.0
    Np, Nd = size(rpos)
    intmat = zeros(T, Np, Np)
    for i in Base.OneTo(Np); for j in Base.OneTo(Np)
        if i == j
            continue
        end
        ri = rpos[i, :]
        rj = rpos[j, :]
        intmat[i,j] = 0.5*Kelbg_repulse(ri, rj; λ11=λ11)
        #ϕ += 0.5*Kelbg_repulse(ri, rj; λ12=λ12)
    end; end
    return sum(intmat)
end


"""
电子和核之间的相互作用
"""
function Kelbg_attract(r1::Vector{T}, r2::Vector{T}; λ22=_WLNeur, Cne=_NeurCoul) where {T}
    return -Cne*_Kelbg_potential(r1-r2, λ22)
    #dis = @. (r1 - r2)^2
    #rabs = sqrt(sum(dis))
    #rrel = rabs/λ22
    #c2 = -exp(-rrel^2)
    #c3 = sqrt(π)*rrel*(1-erf(rrel))
    #if isapprox(0.0, (1 + c2), atol=eps())
    #    return -Cne*sqrt(π)/λ22
    #else
    #    return -Cne*(1 + c2 + c3) / rabs
    #end
end


"""
所有电子和核之间的作用
"""
function Kelbg_attract_pairs(epos::Matrix{T}, npos::Matrix{T}; λ22=_WLNeur, Cne=_NeurCoul) where {T}
    ϕ = 0.0
    Np, Nd = size(epos)
    Nc, Nd = size(npos)
    intmat = zeros(T, Np, Nc)
    for i in Base.OneTo(Np); for j in Base.OneTo(Nc)
        ri = epos[i, :]
        rj = npos[j, :]
        intmat[i,j] = Kelbg_attract(ri, rj; λ22=λ22, Cne=Cne)
        #ϕ += 0.5*Kelbg_repulse(ri, rj; λ12=λ12)
    end; end
    return sum(intmat)
end


"""
所有电子和核之间的作用
"""
function Kelbg_attract_pairs2(epos::Matrix{T}, npos::Matrix{T}; λ22=_WLNeur, Cne=_NeurCoul) where {T}
    ϕ = 0.0
    Np, Nd = size(epos)
    Nc, Nd = size(npos)
    intmat = zeros(T, Np, Nc)
    for i in Base.OneTo(Np); for j in Base.OneTo(Nc)
        ri = epos[i, :]
        rj = npos[j, :]
        intmat[i,j] = Kelbg_attract(ri, rj; λ22=λ22, Cne=Cne)
        #ϕ += 0.5*Kelbg_repulse(ri, rj; λ12=λ12)
    end; end
    return intmat
end
